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1. Introduction and conclusions 



In recent years a new approach to the simulation of statistical systems on the lattice has been 
developed which goes under several names: world line or loop gas formalism, all-order strong 
coupling (or hopping parameter) simulations or 'worm' algorithm methods. The key idea is to 
first reformulate the system as the (complete) sum of its strong coupling graphs. This refers to a 
very simple form of the strong coupling expansion that converges as long as the volume remains 
finite. At large correlation length the graphs that have to be included to obtain precise results are 
however of forbiddingly high order for a systematic expansion. Using algorithms similar to those 
proposed in [[]]] and [^] it has become possible on the other hand to estimate the expansion of many 
observables by a Monte Carlo procedure that samples a representative subset of contributions. An 
important bonus is that this approach has been demonstrated to at least in some cases be free 
of critical slowing down or free of sign problems where this is not so with known methods in 
the conventional formulation. The problem of efficiently producing independent long distance 
correlated field configurations is translated into the need of efficiently passing between the relevant 
large strong coupling graphs. This problem seems to be sufficiently different to make progress in 
cases where for instance cluster algorithms in the conventional setup do not work. 

The topic has been reviewed before [||] at Lattice 2008. It also seems to be closely watched by 
the finite-jit QCD community and typically fills a subsection in their reviews, see [Q] and Sourendu 
Gupta's contribution to this conference. These may be consulted in particular in connection with 
progress on the sign problem. In the present contribution we mainly focus on the important success- 
ful extension to non Abelian spin models of the O(A0 and CP(/V — 1) type. While the method here 
is not confined to two dimensions, most tests are conducted there because of asymptotic freedom 
and the possibility to probe deeply into the continuum limit. We also cover progress on fermions 
which unfortunately at present is confined to two dimensions 1 and allows for instance simulations 
of the Gross Neveu model. Finally we mention here that first steps have been made toward the 
treatment of gauge theories, see [§]. 



2. The idea: Ising model as an example 

A good starting point to explain the strategy is the two point correlation in the Ising model 

w u 2- V I{aW= ± i}e^^ gWgCy) <T(«)g(v) Z 2 (u,v) n 

(O(u)O(v)) = 5-= , . , : = . (2.1) 

2^I {ffW=±1} e^>-M-M Z 



Here u,v are sites on a hypercubic periodic lattice of V sites in arbitrary dimension. In ( |2.1| ) we 
emphasize the view of a correlation as a ratio of two partition functions with and without field 
insertions. 

In any finite volume the expansions 2 of Zq,Z2 in powers of j8 are convergent for all values of 
j3. This includes the vicinity of the critical point and all situations where Monte Carlo simulations 
are performed. It will turn out that in general very high orders in /3 are required to realize this 



1 See however ^ for ideas for an approximate method beyond D = 2. 

2 An expansion in powers of tanh j8 instead of /3 would appear more efficient for the Ising model, but would be less 
easy to generalize below. 
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Figure 1 : A contribution in ' 
boundary conditions 



(left graph) and ^2 (right graph). The boundaries are identified due to torus 



convergence and achieve precision. This is not possible in a systematic truncated expansion as there 
are unmanageably many terms or graphs. As in other physical cases it comes to rescue that not all 
terms are needed. A Monte Carlo procedure will instead sample a sufficiently 'important' subset of 
high order terms. Usual systematic strong coupling expansions are restricted to small correlation 
lengths but, on the other hand, allow to take the thermodynamic limit of quantities like Z2/Z0 term 
by term in the expansion. It is through this step that a finite radius of convergence emerges which 
in many cases (certainly in the Ising model) corresponds to a physical phase transition. 
The expansion is set up by using 



e- 



k=Q 



I ^o{x) k o{y) k (2.2) 



for each neighbor pair on each link / = (xy) introducing independent integers k{l) = 0, . . . ,°° on 
all links. For each configuration k = {&(/)} the spins may now be summed over and the partition 
functions are given as 

Zo= £ P*Mw[k], Z 2 (u,v)= £ p* k ®W[k]. (2.3) 

ge% geS%(«,v) 

In this formula the ^-configurations are viewed as graphs of the type shown in Fig. [j]. There are 
k(l) lines on each link. The left graph is in the set % where the spin summation has enforced 
the constraint that each site must be surrounded by an even number of lines and all different such 
graphs make up The set ^(u, v) is visualized in the right graph in Fig. [j] with the red crosses 
at u, v surrounded by an odd number of lines due to the extra spin insertions. The factor 

"v-Vm (24) 



completes the weight implied by (2.2) 



By differentiating InZo we derive the identity for a link (xy) = /' 

P(a(x)a(y)) = ±- £ P^W[k]k(l') = «*(/')»<>. (2-5) 
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As indicated before a typical graph is thus 0(V) in /3 close to the critical point where the left hand 
side is 0(1). 

A direct simulation of the ensemble Zo in the form ( |2.3| ) was tried a long time ago in [§. The 
authors designed a Monte Carlo algorithm that samples graphs in % (mainly) by local deforma- 
tions over plaquettes. They observed critical slowing down comparable to other standard methods. 
Probably mainly for this reason the approach does not seem to have been pursued much further 
at the time. Also the accessibility of physically interesting observables was not obvious in this 
formulation. The two point function for example could in principle be estimated as a product over 
strings of k(l) but this would probably be inefficient at long distance due to a large variance. 

A breakthrough was achieved much later in [[!]] and The essential idea was a joint simula- 
tion of Zo and Z2 in an ensemble with the partition function 

JT= ^^%[i]=£z 2 ( U) v) (2.6) 

ge%_ u,v 

where the sum over without arguments is over graphs with all possible insertion points 

^2 = U KjV Sf 2 (M,v). (2.7) 

Note that the graphs % contributing to Zo are also included as diagonal contributions 3 with u = v. 
Expectation values are now defined as 

m = ^ZP Z ' mW lk}A[g}- (2-8) 



The identity ( |2.5D now reads 

((<V)) 

and summing over all links I = (xy) we measure the internal energy of the original Ising model. In 
addition it is easy to establish the connection for general correlations 

<ff(xMO)> = (2-10) 

which in particular implies ((S U)V )) = %~ l ^ V~ l with the magnetic susceptibility %. The fraction 
of sampled graphs that belongs to % gets smaller toward the critical point but remains larger than 
one out of V. 

In a very simple but useful generalization we include a nonnegative weight p _1 (w — v), into 
our strong coupling ensemble 

3r= J £p^W[k]p-\u-v) => (a(xM0))=p(x) {i ^-;l ) . (2.11) 

ge^ 2 \\°u,v/l 
We adopt the normalization p (0) = 1 and p must respect the lattice periodicity. The advantages of 



this modification will be discussed in sect. 3.3 



Each graph g £ Sfo appear V times in the sum with all possible u ■■ 
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The essential move in a Monte Carlo simulation of the ^2 ensemble with partition function 3? 
is now the following local update step. One may move u to one of its nearest neighbors by shifting 
it over one of the 2D links attached to it. At the same time the k{l) of that link is changed by 
±1 (adding or removing a line of g). Of course similar moves may be made at v, alternatingly or 
picking randomly one of the two insertion points. These allowed moves staying within ^ may be 
used now as proposals for Metropolis acceptance steps. We do not describe any realization here in 
all details, but refer to the literature. Concrete procedures for the j8 expansion discussed here may 
be found in [^] or Algorithms for the tanh/3 expansion of the Ising model with k(l) G {0, 1} 
are discussed in detail in [^], [|10]. 

In these papers it is numerically demonstrated that strong coupling simulations of the Ising 
model have very much reduced and in many cases completely eliminated critical slowing. We 
conclude that it is advantageous to enlarge the graph space from % to ^2 by allowing defects. This 
is true even if we measure in % only as in ( |2.9| ), but we have seen that the 'intermediate' graphs 
contain even more interesting information. In these proceedings we consider the Ising discussion 
only as a preparation for more elaborate models and hence do not review performance results here 
in more detail. 



3. Nonlinear sigma models 
3.1 0(N) 

In [ |ll| ] a generalization of the above strong coupling graph representation to the O(A0 invariant 
nonlinear sigma model has been given. In this case Z2 generalizes to 



Z 2 (w,v) 



l\jd N s8(s 2 -l) 



(3.1) 



where s(x) is an N component unit vector integrated over the sphere. For N = I the Ising model 
is obviously recovered. To generate an expansion in /3 we again expand the Boltzmann factor on 
each link. To then integrate out the spins for each term in this expansion we need the integral 
over the sphere with an arbitrary monomial in the spin as integrand. This information follows by 
differentiation of the generating function with an ,/V component source j 



d N s8{s 2 -\y s =Y J c[n;N]{j-i) n 

n=0 



with coefficients 



c[n;N] 



r(N/2) 



(3.2) 



(3.3) 



2 2n nT(N/2 + n) 

deriving from the expansion of modified Bessel functions. Working out the combinatorics, i.e. the 
multiplicities of each term, we arrive at 4 



N \g\ 



xp (u — v). 



(3.4) 
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(div k)(x)=4 



Figure 2: Schematic view of a graph for the O(A0 model. The divergence div k means the number of 

lines around a site. 



Several elements in this formula require explanation: 

• The graphs ^2 now differ from those in Fig. [I] only in so far that the even number of lines 
surrounding a site are connected pair-wise except two lines ending at u,v. This is visualized 
for a simple case in Fig. ||. 

• By the pairings (corresponding to contractions of the N spin components) a number of \g\ 
closed loops form in the graph g, each contributing a factor N. 

• The weight W collects lc(l) \ factors and the c[n;N] from the site integrations. 

• The symmetry factor 5? [g] generically equals unity. Only if a graph has extra symmetries 
under the exchange of lines then it equals the order of this symmetry group. This is analogous 
to symmetry factors in Feynman diagrams. More details on this subtlety are found in the 



erratum of [11]. 



The relation between the spin correlation and the graph ensemble is changed only by the spin 
contraction 

(s( X )-s(0))=p(x) {{ ^-; )) . (3.5) 



The simulation of ( |3.4| ) requires in addition to the update steps outlined before reroute moves where 
the local line connectivity is changed. It suffices to go to u (or v) and to randomly divorce one of 
the line-pairs passing through 5 . The single line previously ending at u is remarried to one of the 
divorcees with the other one becoming the new single line. This is again employed as a Metropolis 



proposal that is accepted with a well defined probability dictated by the weights in (3.4). By taking 
into account the asymmetric a priori proposal probabilities, also S^[g] is seen to be implemented 
correctly. 

To implement the simulation just sketched, the graph structure including its connectivity has to 
be mapped on a discrete structure in the Computer. This is possible by a linked list. In the language 



4 The symbols &2,W and later 5? are re-used for the different classes of models that we discuss although their 
precise form becomes context dependent in this way. Factors included in W in the references are sometimes pulled out 
and made explicit in this write-up. 

5 If there are no such lines (as in Fig. no move is made. 
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C the configuration can be coded into pointers where each line element residing on a link gets a 
name and points to its successor and predecessor along its closed loop. In addition there must be 
integer variables describing the geometrical embedding of the graph on the lattice. 

The weight N' 8 * may be implemented either exactly or stochastically. In the original publi- 



cation [11] the exact algorithm (R-algorithm) has been described. In this case the /V-dependent 



weight ( p.4| ) is fed into the Metropolis decision. Then it must be known in the reroute step, if the 
passing through line, that is picked for swapping the connectivity, belongs to the line connecting u 
and v or if it is part of a separate closed loop. To make this nonlocal information available requires 
to sometimes travel around one of the closed loops by following the corresponding pointers. As 
typical loop circumferences grow in the continuum limit, an elementary reroute step costs more 



than 0(1) operations. The numerical observation in [ 1 1 ] can now be summarized as follows. There 



is (practically at least) no slowing down in units of iterations corresponding to order 0{V) ele- 
mentary steps. As they cost however slightly more than 0(V) operations there is a small effective 
critical exponent. It was estimated around z ~ 0.3. This refers to the 0(3) model in D = 2 with 
large volume and correlation lengths ^ = 7, . . . , 65 and to the critical D = 3 model at L = 32, 64. 
The name R-algorithm derives from the fact that here N may be taken also to non-integer values 
by continuing the weights. For integer N one may formulate the I-algorithm where the weight 
is incorporated stochastically. Then each closed loop as well as the line between u and v carries an 
integer degree of freedom i = 1,2, . . . ,N that is independently summed over. In the reroute step, 
only lines with the same i can join. Additional update steps are now needed to move the /-labels. 
A minimal way to do this is to randomly assign a new label to the line between u and v after 0(V) 
elementary steps. In a short test Tomasz Korzec has verified that this form of the I-algorithm shows 
very little critical slowing down in the 0(3) model for D = 2 and % in the range mentioned before. 
A more detailed description of the R- versus I-algorithm together with numerical results is given in 



the paper Q12| ] about the loop formulation of the CP(N — 1) model to which we come in sect. |3_3 



3.2 0(3) model with Nienhuis action 



We now imagine to restrict the graph summation in (2.8) or (3.4) to the subclass of graphs 
which obey the constraint k(l) < 1 on all links, which can be easily implemented in the simula- 
tions. In the Ising model this changes the j3 expansion into the tanh/3 expansion. Thus, if we 
accompany the restriction by this substitution we obtain exactly the same correlation and hence 
u — v distribution as before. In general the reduced set of graphs is equivalent to starting from 



Z 2 (w,v) 



fl Jd N s8(s 2 -\) 



n {l + Hx)- S (y)\ 



S\U) 



■s{v) 



(3.6) 



instead of (pT|). This is clearly a new lattice model - hence we rename j8 — > j8 - and the question 
arises if it belongs to the same universality class despite its strange appearance. Such an action has 
been introduced before [|l3|] and studied in detail by Nienhuis [Qj. An exact solution was obtained 
by summing the strong coupling expansion on honeycomb lattices where the loops now cannot 
intersect. The critical region could be reached for N ^ 2 with /3 ^ 1, i.e. a nonnegative weight 



in (3.6) and universality was supported. In our simulations there is absolutely no sign problem 
when taking j3 > 1 where in the original path integral there seems to be a drastic sign problem. We 
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Figure 3: Continuum extrapolation of a step scaling function based on (3.6) 



thus simulated the 0(3) model with k{l) ^ 1 and quickly found that no criticality was reached with 
[5 ^ 1 . To investigate universality we computed a step scaling function [|l^] for the finite volume 
mass gap extracted from time slice correlations. In Fig. || we see that the data points accurately 
extrapolate to the star which is the exact universal answer known fll6[ ] for this case. These runs 
involve values of jS in the range 1.8 ... 3.1. Our careful conclusion is that at least for this special 
case a universal result is reproduced by ([T(J) at a significant precision. 

3.3 CP(tf-l) 

Another class of nonlinear sigma models that are of physical interest are the CP(jV — 1) sys- 
tems. There the spins label one dimensional subspaces in complex space and may be parameterized 
by <p (x) G C N , 1 (x) | = 1 where <p differing by a phase have to be identified. There are two popular 
lattice actions compatible with this structure. One is the explicit gauge field action 



-S[^,U]=l5Y,[U(x,^(x^(x + fi) + U- 1 (x,^(x + fl)^(x) 

Xfi 



(3.7) 



where nearest neighbors are coupled with a U(l) gauge field U(x,n). It is independently integrated 
over without an action of its own. As it can absorb local phase changes to <p the geometric structure 
of the model is respected. A second option is provided by the quartic action 

-5j0] = 2j3,£|0 t (x)0(x + J u)| 2 (3.8) 

which exhibits local U(l) invariance without extra fields. The standard expectation is that these 
actions are in the same universality class and produce the same continuum quantum field theory. 

A convenient way to probe the model is by correlations of the adjoint local density j a (x) = 
<j)^ (x)A"0 (x), where X a are a basis of hermitian traceless N x N matrices normalized by tr(X a X b ) = 
28 ab , i.e. generalized Gell-Mann matrices. Again we focus on the two point function 

Z 2 (u,v) 



Z 



(3.9) 
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Figure 4: Autocorrelation times in the two dimensional CP(3) model in Monte Carlo time units comparable 
to 'sweeps'. 



By steps that generalize those of the Ising and O(A0 cases and which can be found in detail in [ 12] 



we construct a strong coupling representation (first for (3.7)) 



N \s\ 



xp l (u — v). 



(3.10) 



The complex field variables lead to a modified graph structure 



• Each line on a link carries an orientation (arrow) and they are paired at the sites in a way 
respecting the sense of arrows. 

• There are two lines of opposite orientation running between u and v. 

• On each link there is the same number of arrows in either direction. 

The last constraint is a direct consequence of integrating out the U(l) gauge field. In the exponent 
k{l) is the number of lines per orientation. The weight W is again a local product of explicitly 



known [12] terms and 5?\g\ is the symmetry factor. The connection with the adjoint correlation 
can in this case be written as 



(/(0)/W)=pW 



28qb ((5„-V,*)) 

7Y(7Y+1) «SL V » ' 



(3.11) 



If we repeat the construction starting from the quartic action we arrive at exactly the same graph 
structure but obtain a different expression for W and have to replace j8 2 — > fi q . 



It turns out [12] that this model may be simulated by a procedure very similar to the one 
described for the 0(/V) model. For the I-implementation treating the factor stochastically we 
show measured integrated autocorrelation times in Fig. |j. We consider a series of simulations based 



on (3.7) with D = 2, N = 4 and fixing Lj% « 10. The observable K is defined as the right hand side 
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Figure 5: Effective mass as a function of time slice separation t for L = 780 « 10^ in the CP(3) model. 



of ( j2-9j ) (summed over /) which also here is equivalent to the internal energy. The susceptibility % 
and the mass m are extracted via the correlation ( |3.11 ). 



By computing a step scaling function the universality between the two lattice realizations of 
the CP(N — 1) model with actions (3.7) and (3^) has been confirmed at high precision (see Fig. 4 
in [(T2l])forD = 2,Ar = 3. 



In CP(N — 1) models topology is of special interest. The status here is that the extension of 



( p. 10| ) to include a d term is given in fll2[]. A simulation of this modified system remains to be 
tested. Amplitudes in this case are not strictly positive any more, but it is not known for which 
values of 6 this leads to numerically problematic sign fluctuations. 



It is now time to come back to the usage of the free weight p in our simulations. We read ( |3.11| ) 
as follows: If we are able to guess the behavior of the two point function and use this guess for 
p(x), then, up to known factors, the histogram ((5 M _ VjX )) yields the correction factor that turns our 
guess into the exact answer. For a perfect guess, it would be constant, in other words all possible 
separations u — v would occur with equal frequency. Up to autocorrelation effects we would then 
expect equal relative statistical errors at all separations. One can thus trace the exponential decay 
without degrading signal to noise ratio. That this really works in practice is shown in Fig. ||. Here 
p was set to the scalar propagator with mass 10/L. Similar plots are available for Ising and 0(/V) 
simulations. 
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4. Fermions 



In [ |17[ ] the problem of a 2D Majorana- Wilson fermion in an external scalar field on a two 
dimensional torus was formulated as a loop gas and simulated by cluster methods. Urs Wenger 



first proposed [18] that in this formulation this system could alternatively be simulated with the 
'worm' algorithm [Q] which then triggered the independent study JT^]. In this case the relevant 
partition function with two insertions is 



W«fl(v)> 



n- 



(4.1) 



We here integrate over two Grassmann variables per site, P(y — x) is the Wilson projector 4(1 — 
n nYn) if tne link is y = x + n and t, stands for with charge conjugation ^y^^ 7-1 = — 
An external field is present if m = m(x) is not constant. Due to the projector nature of P and the 
Grassmann nilpotency, the expansions of the link factors have only two terms each (k(l) = 0, 1) 
as for the Nienhuis action. Moreover there can be at most two lines adjacent to a site and as a 
consequence loops and the line between u and v cannot intersect. A typical configuration is shown 
in the left panel of Fig. for a massless free fermion. In this case the correlation is given by 



<&(*)«/,((>)> =p(*; 



((8 XM - v d>(k)M a p(k))} 



(4.2) 



where is a well defined and readily computable sign and M a p(k) is a set of 2 x 2 matrices 
which only depends on the directions in which u and v are approached by their connecting line. 
Although the observable can change sign, this two dimensional fermion has no serious sign prob- 
lem, as neither the denominator nor the numerator gets very small. This is because the bulk of 
closed loops that close without winding around the torus are all positive. This is a specialty of two 
dimensional fermions which here appears as the following feature. For each closed loop there is 
the usual fermionic minus sign, but it is canceled by another sign. The latter arises from the trace 
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of a string of Wilson projectors multiplied along the loops. This minus sign can be understood as 
the sign that arises when a spinor is parallel transported around the loop and is rotated by 2%. For 
subtleties related to (anti)periodic boundary conditions of the finite system the reader is referred to 
[19] and to [|20|]. In this publication it is also shown how N species of this loop gas can be coupled 
to represent and simulate the 0(N) invariant Gross Neveu model in a remarkably efficient way. 

The loop representation of the Majorana fermion was generalized from D = 2 to D = 3 in JI^]. 
Here we still have only two Dirac components and the efficient computability is unchanged. In the 
right panel of Fig. |6] we see a configuration from a simulation of a free fermion at m = 0.65 on 
a 12 2 x 24 lattice where we reproduced the propagator very precisely. Making however the mass 
smaller and/or the system much larger we very abruptly encounter the full-blown sign problem. 
The reason is that in D = 3 the positivity of loops only holds as long as they are planar. If non 
planar loops become abundant, the spin factor assumes all values in Z(8) and does not cancel the 
Fermi minus any longer. Instead of Z(8) there would formally be a U(l) in the continuum, which is 
reduced to Z(8) on the lattice along with rotations being reduced to the hypercubic subgroup. This 
fermionic sign problem is presently unsolved. We consider the simple D = 3 Majorana fermion as 
a good laboratory for further thinking. 



5. Triviality of (p 4 theory 

In four dimensions the number one textbook example for a quantum field theory with a self- 
interacting real scalar field is believed to be trivial, i.e. a free field, once the continuum limit is 
taken. While this is rigorously known to be the case for D > 4 and false for D < 4 in the borderline 
case of D = 4 the belief in triviality rests on numerical demonstrations. As a byproduct of the 
strong coupling reformulation discussed here, we have found a very much improved handle on 
such numerical checks for the Ising limit of (p 4 which is the most interesting parameter range 
for triviality. One of the techniques to obtain rigorous bounds in D > 4 has been developed by 
Aizenman JO]]. He uses nothing but the all-order strong coupling form that we have developed in 
sect. ||, called random current representation by him. Quantized currents k(l) flow through the 
links and are conserved mod 2 at sites, with two sources at u and v. By borrowing his replica trick 
and a graph theoretical proposition we could establish ^ the following identity for the Ising model 
at arbitrary D and volume L D 

g R = -^{m R f = 2z D ((X)) { ^ )Mr (5.1) 

In this formula for the usual renormalized coupling g R 

• x is the 2-point susceptibility, 

• %A is the 4-point (connected, symmetric phase) susceptibility, 

• niR is a renormalized mass, for example using the second moment definition, 

• any fixed values of z = m^L defines a renormalization scheme, 

• the simulation samples two independent replica of graphs g,g' with corresponding 

k,k' ,u,u' ,v, v', 
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Figure 7: Cutoff dependence of the renormalized coupling for D = 4, z = 4. By starting integrations of the 
Callan Symanzik equation at the leftmost data point the lines are produced. 

• SC € {0, 1} is an observable computed as follows. SC = 1 holds iff all four defects are in 
one cluster of an auxiliary bond percolation problem. The bond variables in this problem 
are 'off' only on links where k(l) = = k'{l) holds and 'on' on all others. From cluster 
simulations we know how to efficiently compute 3£ . 

Note that ( |Q| ) implies x 4 ^ 0, i.e. the Lebowitz inequality is manifest in this estimator. The 
advantage of our (Aizenman's) method lies in not having to perform a numerical cancellation to 
compute % 4 which avoids a large significance loss. Triviality now amounts to the question whether 
or not gR \ as L/a — » °°. Here for each L/a, j8 is determined by tuning z to the chosen value. 

In [||] a study was made for a relatively small volume z = 2. Although the results for D = 3 , 4, 5 
are consistent with the triviality expectations, in D = 4 there remained some tension in matching 
with the perturbative coupling evolution close to the continuum limit. Probably this must be at- 
tributed to the weakly damped fluctuations of the constant mode. We therefore made another study 
with z = 4 shown in Fig. ^. The data points with L/a = 8, ... ,64 have errors of only about the 
symbol size in spite of only modest CPU time invested on some PCs. The coupling g may be 
identified with gR here. We see a convincing 'convergence' of the perturbative evolution toward 
the data points at small a/L (left side in the plot). If we are willing to conclude that this agreement 
persists for yet smaller a/L then triviality is established for <p 4 (in the Ising limit). 
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